Exploratory Subgroup Identification with ForestSearch

German Breast Cancer Study Group (GBSG) Analysis

Author

Larry F. León

Published

January 30, 2026

1 Introduction

This vignette demonstrates the ForestSearch methodology for exploratory subgroup identification in survival analysis, as described in León et al. (2024) Statistics in Medicine.

1.1 Motivation

In clinical trials, particularly oncology, subgroup analyses are essential for:

  • Evaluating treatment effect consistency across patient populations
  • Identifying subgroups where treatment may be detrimental (harm)
  • Characterizing subgroups with enhanced benefit
  • Informing regulatory decisions and clinical practice

While prespecified subgroups provide stronger evidence, important subgroups based on patient characteristics may not be anticipated. ForestSearch provides a principled approach to exploratory subgroup identification with proper statistical inference.

1.2 Methodology Overview

ForestSearch identifies subgroups through:

  1. Candidate factor selection: Using LASSO and/or Generalized Random Forests (GRF)
  2. Exhaustive subgroup search: Evaluating all combinations up to maxk factors
  3. Consistency-based selection: Applying splitting consistency criteria
  4. Bootstrap bias correction: Adjusting for selection-induced optimism
  5. Cross-validation: Assessing algorithm stability

The key innovation is the splitting consistency criterion: a subgroup is considered “consistent with harm” if, when randomly split 50/50 many times, both halves consistently show hazard ratios ≥ 1.0 (for example if 1.0 represents a meaningful “harm threshold”).

2 Setup

2.1 Load Required Packages

Show code 
library(forestsearch)
library(survival)
library(data.table)
library(ggplot2)
library(gt)
library(grf)
library(policytree)
library(doFuture)
library(doRNG)

# Optional packages for enhanced output
library(patchwork)
library(weightedsurv)

# Set ggplot theme
theme_set(theme_minimal(base_size = 12))

3 Data: German Breast Cancer Study Group Trial

3.1 Study Background

The GBSG trial evaluated hormonal treatment (tamoxifen) versus chemotherapy in node-positive breast cancer patients. Key characteristics:

  • Sample size: N = 686
  • Outcome: Recurrence-free survival time
  • Censoring rate: ~56%
  • Treatment: Hormonal therapy (tamoxifen) vs. chemotherapy

3.2 Data Preparation

Show code 
# Load GBSG data (included in forestsearch package)
df.analysis <- gbsg

# Prepare analysis variables
df.analysis <- within(df.analysis, {
  id <- seq_len(nrow(df.analysis))
  time_months <- rfstime / 30.4375
  grade3 <- ifelse(grade == "3", 1, 0)
  treat <- hormon
})

# Define variable roles
confounders.name <- c("age", "meno", "size", "grade3", "nodes", "pgr", "er")
outcome.name <- "time_months"
event.name <- "status"
id.name <- "id"
treat.name <- "hormon"

# Display data structure
cat("Sample size:", nrow(df.analysis), "\n")
Sample size: 686
Show code 
cat("Events:", sum(df.analysis[[event.name]]), 
    sprintf("(%.1f%%)\n", 100 * mean(df.analysis[[event.name]])))
Events: 299 (43.6%)
Show code 
cat("Baseline factors:", paste(confounders.name, collapse = ", "), "\n")
Baseline factors: age, meno, size, grade3, nodes, pgr, er

3.3 Baseline Characteristics

Show code 
create_summary_table(
  data = df.analysis,
  treat_var = treat.name,
  table_title = "GBSG Baseline Characteristics by Treatment Arm",
  vars_continuous = c("age", "nodes", "size", "er", "pgr"),
  vars_categorical = c("grade", "meno"),
  font_size = 12
)
GBSG Baseline Characteristics by Treatment Arm
Characteristic Control (n=440) Treatment (n=246) P-value1 SMD2
age Mean (SD) 51.1 (10.0) 56.6 (9.4) <0.001 0.57
nodes Mean (SD) 4.9 (5.6) 5.1 (5.3) 0.665 0.03
size Mean (SD) 29.6 (14.4) 28.8 (14.1) 0.470 0.06
er Mean (SD) 79.7 (124.2) 125.8 (191.1) <0.001 0.30
pgr Mean (SD) 102.0 (170.0) 124.3 (249.7) 0.213 0.11
grade 0.273 0.06
1 48 (10.9%) 33 (13.4%)
2 281 (63.9%) 163 (66.3%)
3 111 (25.2%) 50 (20.3%)
meno <0.001 0.27
0 231 (52.5%) 59 (24.0%)
1 209 (47.5%) 187 (76.0%)
1 P-values: t-test for continuous, chi-square/Fisher's exact for categorical/binary variables
2 SMD = Standardized mean difference (Cohen's d for continuous, Cramer's V for categorical)

3.4 Kaplan-Meier Analysis (ITT Population)

Show code 
# Prepare counting process data for KM plot
dfcount <- df_counting(
  df = df.analysis,
  by.risk = 6,
  tte.name = outcome.name,
  event.name = event.name,
  treat.name = treat.name
)

# Plot with confidence intervals and log-rank test
plot_weighted_km(
  dfcount,
  conf.int = TRUE,
  show.logrank = TRUE,
  ymax = 1.05,
  xmed.fraction = 0.775,
  ymed.offset = 0.125
)

The ITT Cox hazard ratio estimate is approximately 0.69 (95% CI: 0.54, 0.89), suggesting an overall benefit for hormonal therapy.

4 Preliminary Analysis: Generalized Random Forests

Before running ForestSearch, we can use GRF to explore potential treatment effect heterogeneity and identify candidate factors.

Show code 
t0 <- proc.time()

grf_est <- grf.subg.harm.survival(
  data = df.analysis,
  confounders.name = confounders.name,
  outcome.name = outcome.name,
  event.name = event.name,
  id.name = id.name,
  treat.name = treat.name,
  maxdepth = 2,
  n.min = 60,
  dmin.grf = 12,
  frac.tau = 0.6,
  details = TRUE
)
tau, maxdepth = 46.75811 2 
   leaf.node control.mean control.size control.se depth
1          2         6.49        82.00       3.34     1
2          3        -4.10       604.00       1.06     1
11         4        -7.90       112.00       2.81     2
21         5         3.86       177.00       1.87     2
4          7        -5.89       356.00       1.33     2

Selected subgroup:
  leaf.node control.mean control.size control.se depth
1         2         6.49        82.00       3.34     1

GRF subgroup found
Terminating node at max.diff (sg.harm.id):
[1] "er <= 0"

All splits:
[1] "er <= 0"   "age <= 50" "age <= 43"
Show code 
timings$grf <- (proc.time() - t0)["elapsed"]
Show code 
# Display policy trees
# leaf1 = recommend control, leaf2 = recommend treatment
par(mfrow = c(1, 2))
plot(grf_est$tree1, leaf.labels = c("Control", "Treat"), main = "Depth 1")
Show code 
plot(grf_est$tree2, leaf.labels = c("Control", "Treat"), main = "Depth 2")
Show code 
par(mfrow = c(1, 1))

GRF identifies estrogen receptor status (ER) as a key factor, with ER ≤ 0 suggesting potential harm from hormonal therapy.

5 ForestSearch Analysis

5.1 Parallel Processing Configuration

ForestSearch supports parallel processing for computationally intensive operations (bootstrap, cross-validation).

Show code 
# Detect available cores (use fewer for CRAN checks)
n_cores <- min(parallel::detectCores() - 1, 4)
cat("Using", n_cores, "cores for parallel processing\n")
Using 4 cores for parallel processing
Show code 
# Configure parallel backend
registerDoFuture()
registerDoRNG()

5.2 Running ForestSearch

ForestSearch performs an exhaustive search over candidate subgroup combinations with up to maxk factors. Key parameters:

Parameter Value Description
hr.threshold 1.25 Minimum HR for consistency evaluation
hr.consistency 1.0 Minimum consistency rate for candidates
pconsistency.threshold 0.90 Required consistency for selection
maxk 2 Maximum factors in subgroup definition
n.min 60 Minimum subgroup sample size
d0.min, d1.min 12 Minimum events per treatment arm
Show code 
t0 <- proc.time()

fs <- forestsearch(
  df.analysis,
  confounders.name = confounders.name,
  outcome.name = outcome.name,
  treat.name = treat.name,
  event.name = event.name,
  id.name = id.name,
  # Threshold parameters (per León et al. 2024)
  hr.threshold = 1.25,
  hr.consistency = 1.0,
  pconsistency.threshold = 0.90,
  stop_threshold = 0.95,
  # Search configuration
  sg_focus = "hrMaxSG",
  max_subgroups_search = 10,
  use_twostage = TRUE,
  # Factor selection
  use_grf = TRUE,
  use_lasso = TRUE,
  cut_type = "default",
  # Subgroup constraints
  maxk = 2,
  n.min = 60,
  d0.min = 12,
  d1.min = 12,
  # Consistency evaluation
  fs.splits = 400,
  # Parallel processing
  parallel_args = list(
    plan = "multisession",
    workers = n_cores,
    show_message = TRUE
  ),
  # Output options
  showten_subgroups = TRUE,
  details = TRUE,
  plot.sg = TRUE
)
Note: stop_threshold reset to NULL for sg_focus = 'hrMaxSG'.
Early stopping disabled when HR is prioritized in selection;
Also consider increasing max_subgroups_search if search does not appear exhaustive.


=== Two-Stage Consistency Evaluation Enabled ===
Stage 1 screening splits: 30 
Maximum total splits: 400 
Batch size: 20 
================================================

GRF stage for cut selection with dmin, tau = 12 0.6 
tau, maxdepth = 46.75811 2 
   leaf.node control.mean control.size control.se depth
1          2         6.49        82.00       3.34     1
2          3        -4.10       604.00       1.06     1
11         4        -7.90       112.00       2.81     2
21         5         3.86       177.00       1.87     2
4          7        -5.89       356.00       1.33     2

Selected subgroup:
  leaf.node control.mean control.size control.se depth
1         2         6.49        82.00       3.34     1

GRF subgroup found
Terminating node at max.diff (sg.harm.id):
[1] "er <= 0"

All splits:
[1] "er <= 0"   "age <= 50" "age <= 43"
GRF cuts identified: 3 
  Cuts: er <= 0, age <= 50, age <= 43 
# of continuous/categorical characteristics 5 2 
Continuous characteristics: age size nodes pgr er 
Categorical characteristics: meno grade3 
## Prior to lasso: age size nodes pgr er 
#### Lasso selection results 
7 x 1 sparse Matrix of class "dgCMatrix"
                 s0
age     .          
meno    .          
size    0.005433435
grade3  0.178139021
nodes   0.049670523
pgr    -0.001812895
er      .          
Cox-LASSO selected: size grade3 nodes pgr 
Cox-LASSO not selected: age meno er 
### End Lasso selection 
## After lasso: size nodes pgr 
Default cuts included from Lasso: size <= mean(size) size <= median(size) size <= qlow(size) size <= qhigh(size) nodes <= mean(nodes) nodes <= median(nodes) nodes <= qlow(nodes) nodes <= qhigh(nodes) pgr <= mean(pgr) pgr <= median(pgr) pgr <= qlow(pgr) pgr <= qhigh(pgr) 
Categorical after Lasso: grade3 
Factors per GRF: er <= 0 age <= 50 age <= 43 
Initial GRF cuts included er <= 0 age <= 50 age <= 43 
Factors included per GRF (not in lasso) er <= 0 age <= 50 age <= 43 

===== CONSOLIDATED CUT EVALUATION (IMPROVED) =====
Evaluating 16 cut expressions once and caching...
Cut evaluation summary:
  Total cuts:  16 
  Valid cuts:  16 
  Errors:  0 
✓ All 16 factors validated as 0/1
===== END CONSOLIDATED CUT EVALUATION =====

# of candidate subgroup factors= 16 
 [1] "er <= 0"      "age <= 50"    "age <= 43"    "size <= 29.3" "size <= 25"  
 [6] "size <= 20"   "size <= 35"   "nodes <= 5"   "nodes <= 3"   "nodes <= 1"  
[11] "nodes <= 7"   "pgr <= 110"   "pgr <= 32.5"  "pgr <= 7"     "pgr <= 131.8"
[16] "grade3"      
Number of possible configurations (<= maxk): maxk = 2 , # combinations = 528 
Events criteria: control >= 12 , treatment >= 12 
Subgroup search completed in 0.01 minutes
Found 12 subgroup candidate(s)
# of candidate subgroups (meeting all criteria) = 12 
Random seed set to: 8316951 
Two-stage parameters:
  n.splits.screen: 30 
  screen.threshold: 0.763 
  batch.size: 20 
  conf.level: 0.95 
# of unique initial candidates: 12 
# Restricting to top stop_Kgroups = 10 
# of candidates to evaluate: 10 

================================================================================ 
TOP 10 CANDIDATE SUBGROUPS FOR CONSISTENCY EVALUATION
Sorted by: hrMaxSG 
================================================================================ 

Rank   HR        N       Events  K    Subgroup Definition
-------------------------------------------------------------------------------- 
1      2.537     61      34      2    er <= 0 & size <= 35
2      2.222     75      41      2    er <= 0 & pgr <= 32.5
3      2.164     68      38      2    er <= 0 & NOT(age <= 43)
4      2.054     61      35      2    er <= 0 & NOT(size <= 20)
5      1.992     64      34      2    er <= 0 & pgr <= 7
6      1.951     82      45      1    er <= 0
7      1.710     72      39      2    grade3 & pgr <= 7
8      1.707     71      36      2    age <= 50 & pgr <= 7
9      1.530     177     55      2    age <= 50 & NOT(age <= 43)
10     1.397     142     72      2    age <= 50 & pgr <= 32.5
-------------------------------------------------------------------------------- 
Parallel config: workers = 4 , batch_size = 8 
Batch 1 / 2 : candidates 1 - 8 
Batch 2 / 2 : candidates 9 - 10 
Evaluated 10 of 10 candidates (complete) 
6 subgroups passed consistency threshold

*** Subgroup found: {er <= 0} 
% consistency criteria met= 0.92 
SG focus = hrMaxSG 
Seconds and minutes forestsearch overall = 4.322 0.072 
Consistency algorithm used: twostage
Show code 
plan("sequential")
timings$forestsearch <- (proc.time() - t0)["elapsed"]

cat("\nForestSearch completed in", 
    round(timings$forestsearch, 1), "seconds\n")

ForestSearch completed in 4.3 seconds

5.3 ForestSearch Results

5.3.1 Identified Subgroups

Show code 
# Generate results tables
res_tabs <- sg_tables(fs, ndecimals = 3, which_df = "est")

# Display top subgroups meeting criteria
res_tabs$sg10_out
Identified Subgroups
Two-factor subgroups (maxk=2)
Factor 1 Factor 2 N Events E1 HR Pcons
{er <= 0} 82 45 16 1.951 0.920
{er <= 0} {pgr <= 32.5} 75 41 16 2.222 0.960
{er <= 0} !{age <= 43} 68 38 14 2.164 0.950
{er <= 0} {pgr <= 7} 64 34 13 1.992 0.920
{er <= 0} {size <= 35} 61 34 15 2.537 0.990
{er <= 0} !{size <= 20} 61 35 12 2.054 0.940
Search Configuration: Single-factor candidates (L) = 32; Maximum combinations evaluated = 528; Search depth (maxk) = 2
Search Results: Candidate subgroups found = 12; Maximum HR estimate = 2.54
Note: E1 = events in treatment arm; Pcons = consistency proportion

5.3.2 Treatment Effect Estimates

Show code 
# ITT and subgroup estimates
res_tabs$tab_estimates
Treatment Effect Estimates
Training data estimates
Subgroup n n1 events m1 m0 RMST HR (95% CI)
ITT 686 (100.0%) 246 (35.9%) 299 (43.6%) 66.3 50.2 7.8 0.69 (0.54, 0.89)
Questionable1 82 (12.0%) 26 (31.7%) 45 (54.9%) 22.9 43.7 -14 1.95 (1.04, 3.67)
Recommend 604 (88.0%) 220 (36.4%) 254 (42.1%) 66.7 52.6 9.3 0.61 (0.47, 0.80)
1 Identified subgroup : {er <= 0}

5.3.3 Identified Subgroup Definition

Show code 
cat("Identified subgroup (H):", paste(fs$sg.harm, collapse = " & "), "\n")
Identified subgroup (H): {er <= 0}
Show code 
cat("Subgroup size:", sum(fs$df.est$treat.recommend == 0), 
    sprintf("(%.1f%% of ITT)\n", 
            100 * mean(fs$df.est$treat.recommend == 0)))
Subgroup size: 82 (12.0% of ITT)

ForestSearch identifies Estrogen ≤ 0 (ER-negative) as the subgroup with potential harm. This is biologically plausible: tamoxifen is a selective estrogen receptor modulator with limited efficacy in ER-negative tumors.

6 Bootstrap Bias Correction

6.1 Rationale

Cox model estimates from identified subgroups are upwardly biased due to the selection process (subgroups are selected because they show extreme effects). Bootstrap bias correction addresses this by:

  1. Resampling with replacement
  2. Re-running the entire ForestSearch algorithm
  3. Computing bias terms from bootstrap vs. observed estimates
  4. Applying infinitesimal jackknife variance estimation

6.2 Running Bootstrap Analysis

Show code 
# Number of bootstrap iterations
# Use 500-2000 for production; reduced here for vignette
NB <- 1000

t0 <- proc.time()

fs_bc <- forestsearch_bootstrap_dofuture(
  fs.est = fs,
  nb_boots = NB,
  show_three = FALSE,
  details = TRUE
)
Ystar matrix generated should be 'boots x N': 1000 x 686

ForestSearch parameters for bootstrap iterations:
  - sg_focus: hrMaxSG 
  - maxk: 2 
  - fs.splits: 400 
  - max_subgroups_search: 10 
  - hr.threshold: 1.25 
  - hr.consistency: 1 
  - pconsistency.threshold: 0.9 
  - n.min: 60 
  - use_twostage: TRUE 
  - use_lasso: TRUE 
  - use_grf: TRUE 
  Bootstrap-specific overrides:
  - grf_res: NULL (forces re-selection)
  - grf_cuts: NULL (forces re-selection)
  - parallel_args: sequential (prevents nested parallelism)
  - details: FALSE (suppressed in workers)
  - plot.sg: FALSE
  - plot.grf: FALSE

=== Bootstrap Analysis Complete ===
Success rate: 86.8% (868/1000)

H (Questionable) Estimates:
  Unadjusted:       1.95 (1.04,3.67) 
  Bias-corrected:  1.52 (0.8,2.91) 

Hc (Recommend) Estimates:
  Unadjusted:       0.61 (0.47,0.8) 
  Bias-corrected:  0.64 (0.45,0.92) 
===================================
Show code 
plan("sequential")
timings$bootstrap <- (proc.time() - t0)["elapsed"]

cat("\nBootstrap completed in", 
    round(timings$bootstrap / 60, 1), "minutes\n")

Bootstrap completed in 13 minutes

6.3 Bootstrap Summary and Diagnostics

Show code 
# Comprehensive summary with diagnostics
summaries <- summarize_bootstrap_results(
  sgharm = fs$sg.harm,
  boot_results = fs_bc,
  create_plots = TRUE,
  est.scale = "hr"
)

===============================================================
           BOOTSTRAP ANALYSIS SUMMARY                          
===============================================================

IDENTIFIED SUBGROUP:
-------------------------------------------------------------
  H: {er <= 0}

BOOTSTRAP SUCCESS METRICS:
-------------------------------------------------------------
  Total iterations:              1000
  Successful subgroup ID:        868 (86.8%)
  Failed to find subgroup:       132 (13.2%)

TIMING ANALYSIS:
-------------------------------------------------------------
Overall:
  Total bootstrap time:          12.95 minutes (0.22 hours)
  Average per iteration:         0.01 min (0.8 sec)
Show code 
# Display bias-corrected estimates table
summaries$table
Treatment Effect by Subgroup
Bootstrap bias-corrected estimates (1000 iterations)
Subgroup
Sample Size
Survival
Treatment Effect
N NT Events MedT MedC RMSTd HR
(95% CI)†1		 HR‡
(95% CI)2
Qstnbl3 82 (12.0%) 26 (31.7%) 45 (54.9%) 22.9 43.7 -14 1.95 (1.04, 3.67) 1.52 (0.8,2.91)
Recmnd 604 (88.0%) 220 (36.4%) 254 (42.1%) 66.7 52.6 9.3 0.61 (0.47, 0.80) 0.64 (0.45,0.92)
1 Unadjusted HR: Standard Cox regression hazard ratio with robust standard errors
2 Bias-corrected HR: Bootstrap-adjusted estimate using infinitesimal jackknife method (1000 iterations). Corrects for optimism in subgroup selection.
3 Identified subgroup: {er <= 0}
Note: Med = Median survival time (months). RMSTd = Restricted mean survival time difference. Subgroup identified in 86.8% of bootstrap samples.

6.4 Kaplan-Meier by Identified Subgroups

Show code 
 km_result <- plot_sg_weighted_km(
   fs.est = fs,
   outcome.name = "time_months",
   event.name = "status",
   treat.name = "hormon",
   show.logrank = FALSE,
   conf.int = TRUE,
   by.risk = 12,
   show.cox = FALSE, show.cox.bc = TRUE,
   fs_bc = fs_bc,
   hr_bc_position = "topright"
 )

Kaplan-Meier survival curves by identified subgroup

Note: Identified subgroup: {er <= 0}. HR(bc) = bootstrap bias-corrected hazard ratio. Medians [95% CI] for arms are un-adjusted.

6.4.1 Event Count Summary

Low event counts can lead to unstable HR estimates. This summary helps identify potential issues:

Show code 
# note that default required minimum events is 12 for subgroup candidate
# Here we evaluate frequency of subgroup candidates in bootstrap samples less than 15
event_summary <- summarize_bootstrap_events(fs_bc, threshold = 15)

=== Bootstrap Event Count Summary ===
Total bootstrap iterations: 1000
Event threshold: <15 events

ORIGINAL Subgroup H on BOOTSTRAP samples:
  Control arm <15 events: 0 (0.0%)
  Treatment arm <15 events: 0 (0.0%)
  Either arm <15 events: 0 (0.0%)

ORIGINAL Subgroup Hc on BOOTSTRAP samples:
  Control arm <15 events: 0 (0.0%)
  Treatment arm <15 events: 0 (0.0%)
  Either arm <15 events: 0 (0.0%)

NEW Subgroups found: 868 (86.8%)

NEW Subgroup H* on ORIGINAL data:
  Control arm <15 events: 21 (2.4% of successful)
  Treatment arm <15 events: 89 (10.3% of successful)
  Either arm <15 events: 100 (11.5% of successful)

NEW Subgroup Hc* on ORIGINAL data:
  Control arm <15 events: 0 (0.0% of successful)
  Treatment arm <15 events: 0 (0.0% of successful)
  Either arm <15 events: 0 (0.0% of successful)

6.4.2 Bootstrap Diagnostics

Show code 
# Quality metrics
summaries$diagnostics_table_gt
Bootstrap Diagnostics Summary
Analysis of 1000 bootstrap iterations
Category Metric Value
Success Rate Total iterations 1000
Successful 868 (86.8%)
Failed 132 (13.2%)
Success rating Good
Subgroup H (Questionable) Observed HR 1.951
Bias-corrected HR 1.521
Bootstrap CV (%) 101.5%
N estimates 868
Subgroup Hc (Recommend) Observed HR 0.615
Bias-corrected HR 0.644
Bootstrap CV (%) 48.5%
N estimates 868

6.4.3 Subgroup Agreement

How consistently does bootstrap identify the same subgroup?

Show code 
# Agreement with original analysis
if (!is.null(summaries$subgroup_summary$original_agreement)) {
  summaries$subgroup_summary$original_agreement
}
                            Metric       Value
                            <char>      <char>
1:      Total bootstrap iterations        1000
2:           Successful iterations         868
3: Failed iterations (no subgroup)         132
4:                                            
5:    Original subgroup definition   {er <= 0}
6:       Exact match with original  91 (10.5%)
7:         Different from original 777 (89.5%)
8:   Partial match (shared factor) 109 (12.6%)
Show code 
# Factor presence across bootstrap iterations
if (!is.null(summaries$subgroup_summary$factor_presence)) {
  summaries$subgroup_summary$factor_presence
}
  Rank Factor Count   Percent
6    1    pgr   414 47.695853
1    2    age   374 43.087558
2    3     er   342 39.400922
7    4   size   275 31.682028
5    5  nodes    90 10.368664
3    6 grade3    87 10.023041
4    7   meno    47  5.414747

6.4.4 Bootstrap Distributions

Show code 
if (!is.null(summaries$plots)) {
  summaries$plots$H_distribution + summaries$plots$Hc_distribution
}

7 Cross-Validation

Cross-validation assesses the stability of the ForestSearch algorithm. Two approaches are available:

7.1 K-Fold Cross-Validation

Show code 
# 10-fold CV with multiple iterations
# Use Ksims >= 50 for production
Ksims <- 100

t0 <- proc.time()

fs_kfold <- forestsearch_tenfold(
  fs.est = fs,
  sims = Ksims,
  Kfolds = 10,
  details = FALSE,
  parallel_args = list(
    plan = "multisession",
    workers = n_cores,
    show_message = TRUE
  )
)

plan("sequential")
timings$kfold <- (proc.time() - t0)["elapsed"]

# Summary metrics
print(fs_kfold$find_summary)
       Any      Exact At least 1       Cov1       Cov2  Cov 1 & 2 Cov1 exact 
       0.6        0.4        0.4        0.4         NA        0.0        0.4 
Cov2 exact 
        NA
Show code 
print(fs_kfold$sens_summary)
   sens_H   sens_Hc     ppv_H    ppv_Hc 
0.4451220 0.9536424 0.5636029 0.9259563

7.2 Out-of-Bag (N-Fold) Cross-Validation

N-fold CV (leave-one-out) provides the most rigorous stability assessment:

Show code 
t0 <- proc.time()

fs_OOB <- forestsearch_Kfold(
  fs.est = fs,
  details = TRUE,
  Kfolds = nrow(df.analysis),  # N-fold = leave-one-out
  parallel_args = list(
    plan = "multisession",
    workers = n_cores,
    show_message = TRUE
  )
)
Cross-validation setup:
  - Observations: 686 
  - Folds: 686 
  - Fold sizes (range): 1-1 

ForestSearch parameters for CV folds:
  - sg_focus: hrMaxSG 
  - maxk: 2 
  - fs.splits: 400 
  - max_subgroups_search: 10 
  - hr.threshold: 1.25 
  - hr.consistency: 1 
  - pconsistency.threshold: 0.9 
  - n.min: 60 
  - use_twostage: TRUE 
  - use_lasso: TRUE 
  - use_grf: TRUE 
  - (per-fold parallel: sequential)
  - (per-fold details: FALSE)
  - (per-fold plot.sg: FALSE)

Cross-validation complete:
  - Time: 9.85 minutes
  - Subgroup found in 97.8 % of folds
Any found: 0.9781341 
Exact match: 0.9752187 
At least 1 match: 0.9752187 
Cov 1 any: 0.9752187 
Cov 2 any: NaN 
Cov 1 and 2 any: 0 
Cov 1 exact: 0.9752187 
Cov 2 exact: NaN 
Agreement (sens, ppv) in H and Hc: 0.8292683 1 1 0.9773463
Show code 
plan("sequential")
timings$oob <- (proc.time() - t0)["elapsed"]

# Summarize OOB results


cv_out <- forestsearch_KfoldOut(
  res = fs_OOB,
  details = FALSE,
  outall = TRUE
)

tables <- cv_summary_tables(cv_out)

tables$estimates_table
Cross-Validation Subgroup Estimates
Comparison of Full-Analysis vs K-fold CV Predictions
Population / Subgroup Definition N nE EventsE EventsC RMST Diff HR (95% CI)
Overall (ITT) ITT 686 (100.0%) 246 (35.9%) 66.3 50.2 7.8 0.69 (0.54, 0.89)
Full-Analysis: Questionable (H)1 Not recommend 82 (12.0%) 26 (31.7%) 22.9 43.7 -14 1.95 (1.04, 3.67)
K-fold CV: Questionable (H)2 Not recommend 68 (9.9%) 16 (23.5%) 38.9 42.9 2.9 0.85 (0.37, 1.97)
Full-Analysis: Recommend (Hc) Recommend 604 (88.0%) 220 (36.4%) 66.7 52.6 9.3 0.61 (0.47, 0.80)
K-fold CV: Recommend (Hc) Recommend 618 (90.1%) 230 (37.2%) 66.3 55 7 0.70 (0.54, 0.90)
1 Full-Analysis: Subgroup assignment from full-data model
2 K-fold CV: Subgroup assignment predicted from training folds
nE = sample size experimental arm; EventsE/EventsC = events in experimental/control; RMST = restricted mean survival time difference
Show code 
tables$cv_metrics_table
Cross-Validation Agreement Metrics
Agreement between Full-Analysis and K-fold CV Predictions
Metric1 Description Value2
Sensitivity (H) Proportion of FA H patients also predicted H in CV 82.9%
Sensitivity (Hc) Proportion of FA Hc patients also predicted Hc in CV 100.0%
PPV (H) Proportion of CV H predictions that match FA H 100.0%
PPV (Hc) Proportion of CV Hc predictions that match FA Hc 97.7%
1 H = Harm/Questionable subgroup; Hc = Complement/Recommend subgroup
2 Sensitivity = True positive rate; PPV = Positive predictive value
Show code 
tables$sg_finding_table
Subgroup Finding Summary
How often the identified subgroup was recovered in CV folds
Metric Description % of Folds
Any Any subgroup found in fold 97.8%
Exact Exact match (both factors) 97.5%
At least 1 At least 1 factor matches 97.5%
Cov1 Covariate 1 found (any threshold) 97.5%
Cov2 Covariate 2 found (any threshold) NaN%
Cov 1 & 2 Both covariates found 0.0%
Cov1 exact Covariate 1 exact match 97.5%
Cov2 exact Covariate 2 exact match NaN%

8 Results Visualization

8.1 Forest Plot

The forest plot summarizes treatment effects across the ITT population, reference subgroups, and identified subgroups with cross-validation metrics.

Show code 
my_theme <- create_forest_theme(base_size = 24, 
footnote_fontsize = 10, cv_fontsize = 10,
  )

# Define reference subgroups for comparison
subgroups <- list(
  age_gt65 = list(
    subset_expr = "age > 65",
    name = "Age > 65",
    type = "reference"
  ),
  age_le65 = list(
    subset_expr = "age <= 65",
    name = "Age ≤ 65",
    type = "reference"
  ),
  pgr_positive = list(
    subset_expr = "pgr > 0",
    name = "PgR > 0",
    type = "reference"
  ),
  pgr_negative = list(
    subset_expr = "pgr <= 0",
    name = "PgR ≤ 0",
    type = "reference"
  )
)

# Create forest plot
# Include fs_kfold and fs_OOB if available for CV metrics
result <- plot_subgroup_results_forestplot(
  fs_results = list(
    fs.est = fs,
    fs_bc = fs_bc,
    fs_OOB = fs_OOB,
    fs_kfold = fs_kfold
  ),
  df_analysis = df.analysis,
  subgroup_list = subgroups,
  outcome.name = outcome.name,
  event.name = event.name,
  treat.name = treat.name,
  E.name = "Hormonal",
  C.name = "Chemo",
  ci_column_spaces = 25,
  xlog = TRUE,
  theme = my_theme
)

# Display
#plot(result$plot)

#render_forestplot(result, width = 6, height = 8)

# Option 2: Custom sizing
render_forestplot(result)

Subgroup forest plot including identified subgroups
Show code 
# grid::grid.newpage()
# grid::pushViewport(grid::viewport(width = 0.95, height = 0.95))
# grid::grid.draw(result$plot)
# grid::popViewport()


# library(gridExtra)
# grid::grid.newpage()
# grid.arrange(result$plot, newpage = FALSE)

8.1.1 KM Difference plots: ITT and subgroups

The solid black line denotes the ITT Kaplan-Meier treatment difference estimates along with 95%95% CIs (the grey shaded region). K-M differences corresponding to subgroups are displayed.

Show code 
# # Core display of ITT and identified subgroups
# plot_km_band_forestsearch(
#  df = df.analysis,
#    fs.est = fs,
#  outcome.name = outcome.name,
#    event.name = event.name,
#    treat.name = treat.name
# )

# Add additional subgroups along with ITT and identified subgroups
ref_sgs <- list(
age_young = list(subset_expr = "age < 65", color = "brown"),
age_old = list(subset_expr = "age >= 65", color = "orange")
)

plot_km_band_forestsearch(
 df = df.analysis,
   fs.est = fs,
 ref_subgroups = ref_sgs,
 outcome.name = outcome.name,
   event.name = event.name,
   treat.name = treat.name,
 draws_band = 1000
)

Show code 
# # Example with more subgroups
# ref_sgs <- list(
# pgr_positive = list(subset_expr = "pgr > 0", color ="green"),
# pgr_negative = list(subset_expr = "pgr <= 0", color = "purple"),
# age_young = list(subset_expr = "age < 65", color = "brown"),
# age_old = list(subset_expr = "age >= 65", color = "orange")
# )
# plot_km_band_forestsearch(
#  df = df.analysis,
#    fs.est = fs,
#  ref_subgroups = ref_sgs,
#  outcome.name = outcome.name,
#    event.name = event.name,
#    treat.name = treat.name
# )

9 Summary and Interpretation

9.1 Key Findings

Show code 
# Extract key results
cat("=" %>% rep(60) %>% paste(collapse = ""), "\n")
============================================================
Show code 
cat("FORESTSEARCH ANALYSIS SUMMARY\n")
FORESTSEARCH ANALYSIS SUMMARY
Show code 
cat("=" %>% rep(60) %>% paste(collapse = ""), "\n\n")
============================================================
Show code 
cat("Dataset: GBSG (N =", nrow(df.analysis), ")\n")
Dataset: GBSG (N = 686 )
Show code 
cat("Outcome: Recurrence-free survival\n\n")
Outcome: Recurrence-free survival
Show code 
cat("ITT Analysis:\n")
ITT Analysis:
Show code 
cat("  HR (95% CI): 0.69 (0.54, 0.89)\n\n")
  HR (95% CI): 0.69 (0.54, 0.89)
Show code 
cat("Identified Subgroup (H):\n")
Identified Subgroup (H):
Show code 
cat("  Definition:", paste(fs$sg.harm, collapse = " & "), "\n")
  Definition: {er <= 0}
Show code 
cat("  Size:", sum(fs$df.est$treat.recommend == 0), 
    sprintf("(%.1f%%)\n", 100 * mean(fs$df.est$treat.recommend == 0)))
  Size: 82 (12.0%)
Show code 
cat("  Unadjusted HR:", sprintf("%.2f", exp(fs$grp.consistency$out_sg$result$hr[1])), "\n")
  Unadjusted HR: 7.04
Show code 
cat("\nComplement Subgroup (Hc):\n")

Complement Subgroup (Hc):
Show code 
cat("  Size:", sum(fs$df.est$treat.recommend == 1),
    sprintf("(%.1f%%)\n", 100 * mean(fs$df.est$treat.recommend == 1)))
  Size: 604 (88.0%)

9.2 Clinical Interpretation

The ForestSearch analysis identifies estrogen receptor-negative (ER ≤ 0) patients as a subgroup with potential lack of benefit from hormonal therapy.

Biological plausibility: Tamoxifen is a selective estrogen receptor modulator. Its efficacy depends on ER expression. The finding that ER-negative patients may not benefit is consistent with:

  • Mechanistic understanding of tamoxifen action
  • Meta-analyses showing no tamoxifen benefit in ER-negative breast cancer
  • Clinical guidelines recommending tamoxifen primarily for ER-positive tumors

Caveats:

  1. This is an exploratory analysis requiring independent validation
  2. The bias-corrected estimates have wider confidence intervals
  3. Cross-validation metrics should be evaluated for algorithm stability

9.3 Computational Timing

Show code 
timings$total <- (proc.time() - t_vignette_start)["elapsed"]

timing_df <- data.frame(
  Analysis = c("GRF", "ForestSearch", "Bootstrap", "Total"),
  Seconds = c(
    timings$grf,
    timings$forestsearch,
    timings$bootstrap,
    timings$total
  )
)
timing_df$Minutes <- timing_df$Seconds / 60

gt(timing_df) |>
  tab_header(title = "Computational Timing") |>
  fmt_number(columns = c(Seconds, Minutes), decimals = 1) |>
  cols_label(
    Analysis = "Component",
    Seconds = "Time (sec)",
    Minutes = "Time (min)"
  )
Computational Timing
Component Time (sec) Time (min)
GRF 0.2 0.0
ForestSearch 4.3 0.1
Bootstrap 778.7 13.0
Total 1,899.3 31.7

10 References

León LF, Jemielita T, Guo Z, Marceau West R, Anderson KM (2024). “Exploratory subgroup identification in the heterogeneous Cox model: A relatively simple procedure.” Statistics in Medicine. DOI: 10.1002/sim.10163

11 Session Information

Show code 
sessionInfo()
R version 4.5.1 (2025-06-13)
Platform: aarch64-apple-darwin20
Running under: macOS Tahoe 26.2

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/Los_Angeles
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] weightedsurv_0.1.0      patchwork_1.3.2         doRNG_1.8.6.2          
 [4] rngtools_1.5.2          doFuture_1.2.0          future_1.69.0          
 [7] foreach_1.5.2           policytree_1.2.3        grf_2.5.0              
[10] gt_1.3.0                ggplot2_4.0.1           data.table_1.18.0      
[13] survival_3.8-6          forestsearch_0.0.0.9000

loaded via a namespace (and not attached):
 [1] gtable_0.3.6         shape_1.4.6.1        xfun_0.53           
 [4] htmlwidgets_1.6.4    visNetwork_2.1.4     lattice_0.22-7      
 [7] vctrs_0.6.5          tools_4.5.1          generics_0.1.4      
[10] parallel_4.5.1       tibble_3.3.0         pkgconfig_2.0.3     
[13] Matrix_1.7-3         forestploter_1.1.3   RColorBrewer_1.1-3  
[16] S7_0.2.0             lifecycle_1.0.4      compiler_4.5.1      
[19] farver_2.1.2         stringr_1.6.0        codetools_0.2-20    
[22] litedown_0.7         htmltools_0.5.8.1    sass_0.4.10         
[25] yaml_2.3.10          glmnet_4.1-10        pillar_1.11.1       
[28] iterators_1.0.14     parallelly_1.46.1    commonmark_2.0.0    
[31] tidyselect_1.2.1     digest_0.6.37        stringi_1.8.7       
[34] dplyr_1.1.4          listenv_0.9.1        labeling_0.4.3      
[37] splines_4.5.1        fastmap_1.2.0        grid_4.5.1          
[40] cli_3.6.5            magrittr_2.0.4       DiagrammeR_1.0.11   
[43] base64enc_0.1-3      randomForest_4.7-1.2 future.apply_1.20.1 
[46] withr_3.0.2          scales_1.4.0         rmarkdown_2.30      
[49] globals_0.18.0       gridExtra_2.3        progressr_0.18.0    
[52] evaluate_1.0.5       knitr_1.51           markdown_2.0        
[55] rlang_1.1.6          Rcpp_1.1.0           glue_1.8.0          
[58] xml2_1.4.0           rstudioapi_0.17.1    jsonlite_2.0.0      
[61] R6_2.6.1             fs_1.6.6